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Abstract: In previous works we have demonstrated how the energy distribution of 

massless decay products in two body decays can be used to measure the mass of decaying 
particles. In this work we show how such results can be generalized to the case of multi-body 
decays. The key ideas that allow us to deal with multi-body hnal states are an extension 
of our previous results to the case of massive decay products and the factorization of the 
multi-body phase space. The mass measurement strategy that we propose is distinct from 
alternative methods because it does not require an accurate reconstruction of the entire 
event, as it does not involve, for instance, the missing transverse momentum, but rather 
requires measuring only the visible decay products of the decay of interest. To demonstrate 
the general strategy, we study a supersymmetric model wherein pair-produced gluinos each 
decay to a stable neutralino and a bottom quark-antiquark pair via an off-shell bottom 
squark. The combinatorial background stemming from the indistinguishable visible final 
states on both decay sides can be treated by an “event mixing” technique, the performance 
of which is discussed in detail. Taking into account dominant backgrounds, we are able 
to show that the mass of the gluino and, in favorable cases, that of the neutralino can be 
determined by this mass measurement strategy. 
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1 Introduction and general strategy of the mass measurement 

A stable weakly interacting massive particle (WIMP), with a weak-scale mass, is a well- 
motivated candidate for the dark matter (DM) of the universe. The reason for this is 
that it naturally results in the observed relic density after thermal freeze-out of the early 
universe (for review, see for example, Ref. [1]). In fact, this WIMP DM can also have weak- 
scale interactions with various Standard Model (SM) particles.^ This feature leads to the 
exciting possibility that DM arising in this manner could be detected non-gravitationally^. 
In these cases, a symmetry under which the dark matter is charged must be present in 
order to stabilize this DM from decaying to purely SM final states; the SM particles remain 
uncharged, thus preventing decay. 

While there are many possible non-gravitational probes of such DM, we focus here 
on production of DM at hadron colliders, specifically the Large Hadron Collider (LHC). 
Specifically, many models incorporating DM of the aforementioned type contain particles 
that are not only heavier than the DM and charged under the DM stabilization symmetry, 

^This is especially true if the WIMP DM arises as part of an extension to the SM that has been invoked 
to solve the Planck-weak hierarchy problem. 

^cf. thus far, the only evidence for DM is via its gravitational effects. 
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but also that interact via SM gauge bosons. If these heavier particles (dubbed “parent” 
particles) do interact via say, QCD, they could be copiously produced at hadron colliders, 
which would then be followed by their subsequent decay into the concomitant DM and SM 
particles. By design, the DM particle leaves no visible trace in the particle detector, thus 
its presence in an event is typically inferred from the missing transverse momentum (^r); 
which can be interpreted as a loss of specificity in the kinematic information of the event. 

The primary goal of this paper is to devise a strategy for the simultaneous measurement 
of the masses of the parent and the DM particles in the associated processes despite this 
loss of information. This strategy for the mass measurement also has further applications 
beyond the study of DM particles; it can be applied to any case where a new particle 
decays to a semi-invisible final state. In fact, the invisible particle neither has to be a DM 
candidate, nor has to be absolutely stable - only insofar as its time-of-flight out of the 
detector is concerned. However, for notational simplicity we shall still refer to it as “DM”. 

A full reconstruction of such a decay chain is typically not possible, given that it 
contains an invisible particle. On top of this, due to the DM stabilization symmetry the 
parent particles are typically pair-produced, implying that each event comes with two 
invisible particles. The presence of two invisible particles involves an even greater loss 
of kinematic information from each event, and poses a sizable challenge in the associated 
mass measurement. Methods using the Mt 2 variable and its variations [2-7] have been 
proposed as a solution to overcome this challenge. This class of variables is well-known 
for its usefulness both in measuring the masses of particles (e.g. Refs. [8] and [9]) and in 
isolating new physics signals from their backgrounds (e.g. Ref. [10]). Despite their utility, 
these variables have a possible drawback when aiming at a precise mass measurement: 
they all require information about the total missing momentum. Unfortunately, a precise 
measurement of the missing momentum is often difficult, for instance due to the relatively 
poor reconstruction of the jets that are usually a part of the overall event structure. This 
is an unpleasant feature of missing momentum measurements, especially in those cases in 
which many of the jets that are involved in the measurement of the missing momentum 
are actually not involved in the decay process of interest. Said another way, in general, the 
missing transverse momentum is measured as the opposite of the sum of the momentum 
of all the reconstructed objects (leptons, jets, photons,...) in the event, which means that 
the measurement of the missing momentum is an inherently “global” measurement of said 
event. 

In light of this, we have recently proposed complementary methods for mass measure¬ 
ments which instead use only the energy of the visible particles. The reason to pursue this 
strategy is, of course, that it relies intrinsically on more “local” information, ideally using 
only a subset of the particles coming from a given decay chain. ^ The main idea behind 
the method that we propose is to use the energy spectra of the visible particles. The basic 
result upon which the method is predicated was shown in Ref. [19]: namely, for a massZess 
child from the two-body decay of an unpolarized parent, the peak in the energy spectrum of 

®See also Refs. [11-16] for other recent methods of mass measurement that do not use the missing 
transverse momentum and Ref. [17, 18] for a general review of mass measurement methods. 
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Figure 1. The three-body decay of interest (left panel) and the effective two-body decay (right 
panel) with the mass of the visible system being mab- (cib) denotes the effective visible pseudo¬ 
particle formed by the two visible particles a and b. 

the child (henceforth denoted as “energy-peak”) seen in the laboratory frame (henceforth 
denoted by “laboratory-frame” energy) is the same as the fixed value of its energy in the 
rest frame of the associated parent (henceforth denoted by “rest-frame” energy). The latter 
value is given by a simple relation in terms of the masses of the two massive particles (the 
parent and the other child) involved in the decay, and hence can give information about 
these masses. In a subsequent paper [20], we then applied this observation to measuring the 
unknown masses in the semi-invisible decay of a heavy new particle involving a multi-step 
cascade of two-body decays.^ 

In this paper, we consider instead a smgte-step, ttiree-body, semi-invisible decay of a 
heavier new particle, which we denote as 

B —)■ Aab (1-1) 

where a, b are visible SM particles and A, B are massive new particles with A assumed 
invisible. In order to deal with this specific decay topology, we need to extend the result 
of Ref. [19] to multi-body decays. The key idea is to map a multi-body final state into a 
two-body one, by the factorization of phase-space. In carrying out this mapping, we will 
take particular care in correctly partitioning and grouping the multi-body decay products 
and selecting an appropriate region of the multi-body phase space upon which to apply 
the above two-body result. 

To reduce the multi-body final state to a two body one, we first form a compound 
system made of all of the visible particles, labeled as a and h in eq. (1.1). We denote 
this compound system as {ah) and graphically represent the corresponding partitioning in 
Figure 1. After this partitioning the decay does not yet look like a truly two-body decay 
because the compound system {ah) does not have a fixed mass. The combination of {ah) 
will have its own phase space in invariant mass. This is apparent from the well-known [22] 
recursive formula for the multi-particle phase-space of N particles of masses mi, ...,m 7 v, 
which can be thought of as the sum of many two-body phase-spaces, a single particle on 
its own and the remaining {N — 1) particles clustered into a single object whose mass now 
depends on the momenta and the angles between the N — 1 clustered particles: 

^We also showed that our energy-peak result of Ref. [19] can be used for “counting”’ DM particles in 
decays [21], which is a powerful probe of the DM stabilization symmetry. 
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(/)Ar(mi, ...,mN) = j d^id<j)2 (mAr,/u(mi, ...,mAr_i) • #Ar_i(mi, ...,mN-i)) ■ (1.2) 

Considering each value of the masses that the compound system [ah) can take separately, 
we can regard the N body final state as a weighted sum of a collection of two-body systems, 
each of which is characterized by the mass of the compound system denoted by n and its 
probability d4>]sf-i. This probability, together with the actual squared matrix element of 
the decay, would give the rate of decay in that particular kinematic configuration. In the 
following, we do not assume any knowledge of the matrix element of the decay and we shall 
make no use of these rates; all we will need for our strategy to work is the ability to represent 
the multi-body final state as the sum of the collection of all possible two-body final states. 
The fact that we do not need to know the rate for each possible kinematic configuration 
of the multi-particle final states is a remarkable point of strength of our method; it is 
especially powerful when applied to newly discovered particles, as their matrix elements 
are a priori essentially unknown. 

For the case of a three-body decay, the above outlined procedure gives 

B —)• Aab = {B ^ A {ab)m^f,) , 

rriab 


where the equality should be taken in the sense of an equivalence. We also remark that for 
practical reasons the integral for the phase-space factorization formula has been discretized. 
In this way, we can form a finite number of compound systems {ab)mab of mass mafe±<5 with 
6 <C rriab- This procedure ensures that each of the compound systems has an approximately 
fixed invariant mass, and we can think of it as a “pseudo-particle” having a width of order 
6. This means that we partition, or “slice”, the data according to the total invariant mass 
of a compound particle formed from the a and b particles and apply the result for two-body 
decays to each mass partition as the overall system has been reduced to an effective two- 
body decay. In the rest frame of the parent particle, each partition of the pseudo-particle 
(ab) has energy: 


E 


(ab) 


ml-m\ + ml^ 
2m B 


(1.3) 


based simply on two-body kinematics for decay of B into A and (ab). Using the appro¬ 
priate extension of the result in Ref. [19] to the case mab 7 ^ 0 (see more on this point in 
Section 2), we are able to extract E'^ab) from the laboratory-frame energy distribution of 
that particular (ab) compound particle. We then repeat this procedure for each of the mass 
partitions in the overall range of mab- When plotted versus the fitted data for 
extracted from the energy distributions should lie along a straight line as per eq. (1.3). It 
is straightforward to see that ms can be determined from the slope of this line and that 
mA can be determined from the intercept on the vertical axis once ms has been deter- 
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mined. The available information can be fully utilized in constructing this straight line by 
analyzing the data in all slices.^ 

We remark that, although the characteristic signature of the production of an invisible 
particle A is missing momentum, our method does not make explicit use this quantity and 
yet still offers a way to obtain a measurement of mass of this particle! In other words, any 
specific property of the invisible particle is almost completely irrelevant to our method, 
except for the assumption that there is at least one invisible particle per decay chain. We 
again emphasize that this achievement is remarkable, especially in comparison with other 
mass measurement methods involving such as M'r 2 and its variants. 

Despite the simplicity of the general idea, there are still some potential issues that 
should be properly dealt with in order to successfully execute the strategy outlined above. 
First, the compound particle (ab), the visible child particle from the effective two-body 
decay, has a non-negligible mass; thus, it is essential to generalize the result in Ref. [19] on 
the energy-peak to the case for a massive visible child particle. Our companion paper [30] 
will be devoted to studying how to deal with these massive child particles in more detail. 
In the present work, we shall merely report the final result of this companion paper and 
use this result for our present investigation. Nevertheless, the discussion presented here is 
largely independent of the derivation of this result. 

As mentioned earlier, the DM model under consideration has the parent particles being 
produced in pairs. If both of them decay to the same final state, a combinatorial ambiguity 
arises in attempting to correctly partition and group those particles originating from the 
same parent; multiple pairs can be formed from the final state as seen in the detector, but it 
is not known a priori which is the correct pairing - that is, that the particles in the pairing 
originated from the same decay. This partitioning is a crucial necessity in forming the 
{ab) compound system that plays the role of the child pseudo-particle. For this reason, we 
allot Section 3 to the thorough discussion of the treatment of this combinatorial ambiguity. 
In particular we propose an “event mixing” technique [23] and Refs. [25-29] as a way to 
remove the combinatorial background. 

To illustrate our method we discuss in detail its application to a specific process. As 
a concrete example, we choose the pair production of gluinos in an R-parity conserving 
supersymmetric model. Here the gluinos are assumed to subsequently decay into bb and 
an invisible light stable neutralino via an Ojff-shell bottom squark: 

pp ^ gg ^ bbbbxx- ( 1 - 4 ) 

This scenario is chosen primarly because it has been thoroughly studied in the literature, 
and thus should be familiar and interesting to a large audience. Indeed, this process has 
also been investigated at the LHC [10, 31-35]. In order to provide a fully realistic example, 
and to demonstrate some of the issues that arise in using our method, we shall incorporate 
the relevant Standard Model backgrounds in our analysis as well. We take particular 

®In practice, one could end up not using some of the slices if treating them becomes too problematic, 
e.g. because of backgrounds or sensitivity to the cuts. The fraction of unused slices will in any case be kept 
to a minimum. 
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care in devising cuts for background rejection so that these selections do not affect the 
shape of the energy spectrum near the peak, which is the critical region of interest for 
our energy peak method. Obviously, the optimization of these cuts is a process-dependent 
issue, and so must be evaluated on a case-by-case basis. The goal of our discussion is 
to present potential systematic uncertainties and biases arising from the specific details 
of our method, such as those induced by phase space slicing, event mixing, imperfect 
knowledge of the background, and overly restrictive event selection criteria. We also present 
other complementary observables that enhance the findings obtained using the energy peak 
method, one of which is the kinematic endpoint of the di-jet invariant mass distribution. 

The rest of the paper is organized as follows. We continue in the next section with a 
discussion of a template function used to describe the energy spectrum of a massive child 
particle. Section 3 is devoted to dealing with the combinatorial ambiguity inherent in our 
chosen decay topology. Then in Section 4, we detail our selected signal process and the 
relevant backgrounds, both those from SM processes and the “event mixing” scheme for 
signal combinatorics. Section 5 contains the main results for the mass measurement of 
the aforementioned example process together with a discussion of several opportunities for 
improvement to the method. In Section 6 we present our conclusions and outlook. 

2 A template for the energy spectrum of a massive child particle 

As outlined in the previous section, the essence of our mass measurement technique is to 
fit the data to get the value of E* for each of the fixed masses of the compound system 
{ab) and fit them onto the straight line in eq. (1.3). The mass of the compound system 
{ah), being a system of two particles, is not fixed and in general spans a range fixed by 
the masses of all the particles in the decay. Since we are a priori unaware of the masses of 
the parent and invisible child particles, it is not possible to know whether or not a given 
value of iTiab is small enough in comparison to those unknown masses to justifiably trust 
the validity of our previous results for effectively massless child particles [19], and thus we 
are motivated to extend the finding to massive child particles. 

The primary difficulty in generalizing to a massive child is the potential loss of correla¬ 
tion between the peak location in the laboratory-frame energy distribution of the massive 
child and its energy in the rest frame of the parent. This can readily be seen when consider¬ 
ing the decay of a massive particle B ^ X tf aX the kinematic end point of the phase-space, 
i.e., ms = rnx + For any value of mx, if will be at rest in the B rest frame, hence 
Efp = m^. If each particular event is boosted to the laboratory frame, the energy of 
becomes simply with 7 ^ being the boost of particle B relative to the laboratory. 

This direct linear relationship between E^ and 7 ^ implies that the shape of the energy 
distribution of particle if in the laboratory frame should simply be that of the boost dis¬ 
tribution of particle B. In this case, it is clear that the peak of the energy distribution of 
the massive child if carries essentially no information about the masses; rather, it carries 
information on the most probable boost of particle B. This is contrast to the “invariance” 
that holds for a massless child: the energy-peak in the laboratory frame is the same as the 
rest-frame energy value irrespective of the details in the boost distribution of particle B. 
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For a more formal understanding of this problem, it is instructive to analyze the Lorentz 
transformation of a massive child particle from the rest frame of its parent particle, where 
it has energy-momentum {E*,p*), to the laboratory frame. Given the boost factor 7 of 
the parent particle and the emission angle of the child 9* relative to the boost direction, 
we find the energy of the child particle in the laboratory frame (denoted hy E) to be 

E = E*y{l +P*Pcos9*) , (2.1) 

where we have used p* = I3*E*. We observe that the laboratory-frame energy E becomes 
equal to the rest-frame energy E* only if 

Denoting this 9* as the “reference” angle, we see that any value of cos 9* smaller (larger) 
than eq. (2.2) gives rise to a laboratory-frame energy value E smaller (larger) than E*. 

To ensure the existence of the reference angle, the boost factor should satisfy the 
following relation: 


1 

When this condition is satished, the energy distribution in the laboratory frame is non-zero 
at E = E*, which is a necessary, but not sufficient, condition to have a maximum at E*. 
Obviously, if for some 7 the condition set by eq.(2.3) is not satisfied, then E > E* for all 
cos 9*, which potentially invalidates the statement that the peak in the laboratory-frame 
energy distribution appears at E*. Actually, for the typical boost distributions of parent 
particles produced at hadron colliders, one can see that if any of the boosts of the parent 
particle(s) lie outside of the range given by (2.3), it is then guaranteed that the peak of 
the energy distribution in the laboratory frame will not be located at the rest-frame energy 
value A more rigorous method of characterizing the energy distribution of a massive 
child will be presented in a separate work [30]. We hope that the argument above, while 
not as rigorous as our companion paper, will convince the reader that the maximum boost 
of the parent particle is the key parameter that controls the position of the peak in the 
laboratory-frame energy distribution of a massive child particle. 

On top of affecting the peak position, the overall shape of the energy distribution for 
the massive child is expected to differ from that for the massless child. This means that 
the function used to fit the massless child energy spectra in previous works cannot be used 
in the present work. In order to obtain a suitable description of the massive child energy 
spectrum, we revisit the corresponding discussion for the case with a massless child particle. 
The value of the energy distribution at a given laboratory-frame energy E is given by a 
Lebesque-type integral within the range of 7 values contributing to the E together with 

®The displacement of the maximum with respect to E* may still be small, but strictly speaking, the 
“invariance” that we demonstrated in Ref. [19] for the massless particle will be broken. 
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the associated weight for the 7 [19]. For the case at hand, this range is found by solving 
eq. ( 2 . 1 ) for 7 with cos 8* = ± 1 , which gives: 


7+{E) = 7 *' 



1-{E) = 7*2 


E 

E* 



(2.4) 

(2.5) 


We see that for a massless child particle, i.e., 7 * —>■ 00 , ^+{E) diverges, whereas in the 
same limit 'j-{E) converges to a finite value: 

lim = + ( 2 . 6 ) 

7*->'00 h* E 


In light of eq. (2.6), we can express the energy spectrum for a massless child that we used 
in previous works [19] as 


exp 


(-«^-7i“)(^))- 


(2.7) 


Motivated by the success of this exponential form in the massless case [19, 20] and 
exploiting the identification in the massless limit eq. (2.7), we propose the following ansatz 
on the shape of the laboratory frame energy spectrum of a massive child 


f{E) = N (exp[-u; • 7 -(Fi)] - exp[-u; • 7 +(FI)]) , (2.8) 

where N and w are a normalization factor and the width of the function, respectively. A 
complete evaluation of the accuracy with which this function describes the energy spectrum 
in the laboratory frame of a massive child will be presented in a dedicated work [30]. For 
the purposes of this paper, it will be sufficient to know that this function reproduces our 
ansatz for massless child particles for 7 * —)• 00 . In any case, we will explicitly show that 
this function provides a good description of simulation data for our example process below. 

A comment on the location of the maximum of this function is in order. The maximum 
of this function coincides with E* only in the limit rc —)• 00 , in which the function becomes 
a 6 function. For all finite values of w, the actual location of the maximum of the function 
is slightly larger than E*. However, we empirically observe that for parent particles that 
would typically be produced at colliders, and for 7 * somewhat larger than 1 , the typical 
value of w is large enough that this effect is negligible. Therefore, we expect that eq. (2.8) 
properly describes a large class of energy spectra. Because the peak location and the E* 
parameter are no longer necessarily the same the determinations of the best fit values of 
E* and w are interrelated when fitting the data to the massive template function eq.(2.8). 
The fact that the maximum of the spectrum is a function of both E* and w is an inherent 
feature of our ansatz for the massive child energy spectrum, which did not exist in the 
massless case of Ref. [19] . In this paper we will study the possible effects that arise in our 

^Again, for the case with massless children, E* conforms to the peak irrespective of w. 
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explicit example due to this feature of eq.(2.8). For a fully general investigation of this 
issue we refer to our companion paper [30]. In Sec. 4 we shall fit the energy distribution 
of each mass partition of the pseudo-particle (ab) both with the new template for massive 
children eq. ( 2 . 8 ), and with its massless limit (i.e 7 * —)• 00 ), the latter of which was the 
template employed previously for massless child particles. The comparison of the results 
from these two templates will allow us to demonstrate the necessity of using eq. ( 2 . 8 ) and 
generalizing what we had used in our previous work [19, 20]. 

3 “Event mixing” to estimate the combinatorial background 

As explained in the Introduction and depicted in Figure 1, the success of our strategy 
hinges on correctly identifying the pairs of particles coming from the same decay side. 
If this identification is done correctly, the idea of phase-space factorization can safely be 
applied to reduce the multi-body final state to a two-body final state. Generally speaking, 
the identification of the correct pairs of particles to be grouped together is a tremendously 
difficult task, as we have no systematic way of knowing which particles are related to the 
same decay, and can thus be correctly paired in the analysis. 

One approach to surmounting this challenge is to attempt to correctly identify the 
pairs of particles that come from the same decay by exploiting some preferential kinematic 
correlation between particles that originate from the same decay. This correlation could 
then be translated into a selection criterion that would correctly pair the appropriate 
particles with relatively high accuracy. Of course, this selection process would not be 100% 
effective and would fail to pick the correct pairings for some fraction of events. As a result, 
we would be left with a certain amount of combinatoric pollution from pairings whose 
constituent particles did not come from the same decay. Several event-by-event strategies 
have been developed to identify which pairs of particles come from the same decay (see, 
for example, Refs. [9, 36-39]). It is however generally true that, in order to maximize the 
chance of pairing particles correctly, the kinematic selection criteria which form the basis of 
each method must be rather restrictive, so as to guarantee a sufficient rejection of unwanted 
pairings. Events that pass these highly constrained kinematic criteria will be preferentially 
selected from isolated regions of the kinematic phase space of the scattering, and therefore, 
the kinematic distributions of the final state will be significantly altered by the imposition 
of these criteria. Our method of mass extraction relies critically on the fidelity of the 
energy distribution around the location of the peak; without this, the templates we use to 
extract the masses return biased information. Thus, because the above set of procedures for 
selecting correct decay siblings greatly disturbs the resultant kinematic distributions, they 
are not suitable for use alongside our method. For this reason, we do not even attempt to 
identify the correct pairings in each event; we instead try to obtain the energy distribution 
of the correct pairs without knowing which are the correct pairs event-by-event. 

In order to determine the distribution of a given observable as it would arise from 
picking just the correct pairs, we use an event mixing technique whose basic idea is as 
follows. We consider a scattering 

pp B 1 B 2 (AiOi 6 i)(A 2 a 2 & 2 ), (3.1) 
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where two heavy particles Bi and B 2 are produced and decay to final states of the same 
kind, B —)• Aah, which are labeled to correspond with their respective parent.® In the 
analysis of events of this nature, we follow the procedure laid out in Section 1 for making 
pseudo-particles out of the a and b particles for all possible equivalent pairings (e.g. ai with 
bi, but also oi with 62 and so on.) From these pairs, we obtain a fully inclusive distribution 
of the observable of interest, which in our case here is the total energy of the pair, 

(all pairs). (3.2) 


In order to obtain the distribution stemming only from the pairs of particles coming 
from the same decay, it is sufficient to come up with an estimate of the distribution that 
stems from the pairs that we would like to discard and subtract it from the fully inclusive 
distribution: 

(same-decay pairs) = (all pairs) — (different-decay pairs). (3.3) 


This equality might seem trivial, but is in fact very powerful. To wit, the object that we 
must have for our method to be successful - the distribution of pairings from the same 
decay - which is quite difficult to obtain in itself, is now expressed in terms of two objects 
that are much simpler to obtain. The first piece is obviously attainable, as it is simply 
the distribution from all the pairs that can be formed in each event. The second piece, 
the distribution from pairs not coming from the same decay, can be estimated by the 
“event mixing” technique, which is based on making a distribution from pairs of jets that 
come from different events. The intuition that justifies the usefulness of the event mixing 
technique originates from the fact that pairs ( 0162 ) and ( 02 ^ 1 ) from the same event are 
made of particles which are produced with almost no kinematic correlation. Therefore, it 
seems reasonable to mimic the effect of these “incoherent” pairs with the pairs of particles 
taken from different events, which intuitively have no correlation. More precisely, we can 
see that the phase-space point from which oi and bi originate in the decay of Bi is very 
close to being uncorrelated with the phase-space point from which 02 and 62 originate in the 
decay of i? 2 - This approximately vanishing correlation between the products of different 
decays implies that pairs formed by particle oi taken from one event, and particle 62 taken 
from another event are expected to be statistically equivalent to pairs ( 0162 ), where both 
particles are taken from the same event. This means that the distributions of a quantity 
over a given sample of events obtained from either the (ai 62 )-type incorrect pairings within 
the same event or from pairing ai in one event and 62 in another event are equivalent. This 
implies that we can estimate the distribution stemming from pairs of particles from the 
same decay as 


da 


(same-decay pairs) 


rs»/ 


da 

dEo^h 


(all pairs) 


da 

dEf;^l, 


(different-events pairs). 


(3.4) 


This observation is at the center of the event mixing idea, and we henceforth denote the 
procedure described by the right-hand side of eq. (3.4) as “mixed event subtraction”. 

^Strictly speaking, Ti’s need not be invisible as long as they are distinguishable from Ui and bi, which 
are assumed indistinguishable from each other. 
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The plausibility of the distribution obtained by pairing particles from different decays 
within the same event being equivalent to the distribution arising from pairings between 
different events can be seen intuitively in certain simple cases. To illustrate, imagine having 
an ensemble E of pairs of particles (i^i, B 2 ), 

E = {(SS'), (Sf \ ), \ ),...} , (3.5) 

where the superscripts denote the associated event number. For simplicity, we take particles 
B to be scalars sitting at rest in the laboratory, where they decay B —)■ abA. It is obvious 
that the distributions made from pairs coming from different B particles in the same event 
of the ensemble are the same as the distribution made from pairs coming from B particles 
from different events in the ensemble. In this example, we are simply sampling the phase 
space of the B decay in two different ways, in one case taking kinematic information 
from instances of the decay that happen at the same time and in the other taking that 
information from instances of the decay that are separated in time. 

However, one must note that the situation described above may not correspond to 
reality. As a concrete example, we take the pair of {B^\ B^'^) and assume that the system 
of the two B particles has a center of mass energy 2mB)- In the laboratory frame, 

the phase-space accessible to the decay products of the two B particles depends on y/sl. If 
we consider the jth event, with the intention of mixing particles from the ith and jth events, 
we are forced to confront the fact that the center of mass energy in the jth event, and 
thus the phase space accessible to the jth event’s particles, is different than that of ith event. 
The mismatch in the phase-space accessible to particles in events at different is clearly a 
potential source of error in the identification of the “different-decay” and “different-event” 
distributions, which naturally poses an a threat to the successful application of the event 
mixing technique. 

It is quite difficult to estimate the size of this inaccuracy, though we expect that for 
typical situations at hadron colliders the event mixing technique works quite well. One 
reason for this is because of the small variance of the boosts of the B particles produced in 
typical collisions at hadron colliders. In addition to this effect, however, there may be other 
potential sources of error in the event mixing technique, and a case-by-case study is needed 
to check the performance of the method. Because of this, we take a pragmatic approach in 
the following and apply the event mixing technique to our example while explicitly checking 
the performance of this method for our example process. 

4 Application to the gluino decay 

We now demonstrate how the general strategy detailed above is realized by taking as an 
example a particular gluino decay channel. We first illustrate the signal process and its 
particular characteristics, and then move to discussing the possible backgrounds of this 
signal process. The discussion of these backgrounds is separated into two categories; 1) 
the real background from SM processes, and 2) the systematic background from incorrect 
pairings of the final state particles used in forming the invariant mass and energy distribu¬ 
tions. As we detail the various processes involved, we shall utilize Monte Carlo simulation 
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to generate the relevant event samples, construct and analyze the appropriate kinematic 
data from these samples, and end with a discussion of the effectiveness of our technique 
based on the results of this analysis. 

4.1 Signal process: gluino decay 

We apply the general idea developed in the previous sections for the case of pair-produced 
SUSY gluinos and their subsequent decay into two bottom quarks and a neutralino via a 
three-body decay: 


pp^ gg ^ bbbb + xx (4.1) 

at the 14 TeV LHC. In terms of the notation used in Sectionl, the gluino and the neu¬ 
tralino correspond to particles B and A, respectively, and two visible particles a and b 
are the bottom quark and anti-quark in a decay chain. In reality, the particle detector 
cannot reliably discriminate between bottom and anti-bottom, thus particles a and b in 
this example are considered indistinguishable. 

Though we are using the specific decay above as a concrete example, we emphasize that 
the decay mode and underlying model at hand are chosen only to enable us to demonstrate 
the proposed technique and that the general idea can be applied to multi-body decay 
processes in other models. We also point out that the applicability of the method is not 
affected in any way by the strengthening of bounds on supersymmetric particles, because 
the method is applicable for parent and invisible (child) particles of any mass. To illustrate 
our technique, we choose the masses of these particles to be 

rrig = 1.2 TeV, rriy. = 100 GeV 

with a decoupled bottom squark, and assume that the only decay mode of the gluino is a 
three-body decay in the form of bbx The Monte Carlo signal for our study is simulated 
using MadGraphB vl.4.8 [41] and the structure of the proton is parametrized by the parton 
distribution functions (PDFs) CTEQ611 [42], evaluated with the default renormalization 
and factorization scale settings of MadGraphB. The production cross section of the paired 
gluinos is computed with MadGraphB and is reported in the first column of Table 1. Since 
we assume that all produced gluinos decay into bbx a-s described above, cr{pp —)■ gg) is 
equal to a{pp —t 55 —)• bbbbxx)- 

The neutralinos in the final state of our signal do not interact with the detectors of the 
LHC, resulting in a missing transverse momentum. The four bottom quarks give rise to 
jets of hadrons - particularly, H-hadrons. The particular characterisitcs of the H-hadrons 
in the jets allows us to distinguish this type of jet from other jets that do not originate 
from the bottom quark, and it is possible to see the traces of 6-quark-initiated jets and tag 

®During the completion of this work the limit on the gluino mass given by the LHC experiments has 
risen to about 1.4 TeV for light x [34, 35]. Despite these new limits ruling out the spectrum we consider at 
a 95% confidence level, this spectrum still serves its purpose as an illustration of the technique. It should 
be remarked that we do not expect qualitative differences in the application of our strategy to the mass 
measurement of a heavier but not yet experimentally excluded gluino. 
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46 -|- pT 
(before cuts) 

46 -|- pT 
(after cuts) 

46 -|- Z{^ vv) 
(after cuts) 

ttbb 

(after cuts) 

[fb] 

54.74 

36.53 

0.48 

0.15 


Table 1. The cross sections for the signal process (before and after a set of cuts are imposed) 
and the (main) background processes (only after cuts are imposed). The cuts are described in 
eqs. (4.2)-(4.4), and also include all the identification and the isolation criteria explained in the 
text. The effect of the 6-tagging efficiency is not taken into account by the numbers in this table. 

them in a large fraction of the events. With the requirement that four of the reconstructed 
jets in the final state have this tag, the signal will feature four bottom jets plus missing 
transverse momentum, 46-|-^t- 

Before closing this section, we remark that the chosen example process poses an extra 
challenge in the application of our method. In fact all visible final state particles are 
mdistinguishable, hence there are three different ways to form pseudo-particles from these 
6 -tagged final state particles that must be all be considered in the analysis. Together with 
the SM backgrounds, this can be interpreted as another background, as explained in the 
next section. 

4.2 Backgrounds 

In this section, we discuss the backgrounds relevant to the signal process defined in the 
preceding section. As mentioned earlier, there are two types of backgrounds: those coming 
from Standard Model processes which give rise to the same signature as our supersymmetric 
process, and that which comes from taking the wrong pairings of two visible particles when 
we evaluate both the energy sum and the pair-wise invariant mass. We start by discussing 
the “real” background from the Standard Model and then we discuss the combinatorial 
background. 

4.2.1 Standard Model backgrounds and event selection 

For our collider signature Ah + pT, the following two processes in the SM are identified as 
the major backgrounds: 

pp — 7 - bbbb + Z ^ bbhb -|- l'u and pp —?■ ttbb. 

The Monte Carlo generation of background events is done using the same event generator 
and input PDFs as those for the signal events. Since the detector signature from these 
interactions is exactly the same as the one used in our earlier work Ref. [20], we adopt 
a similar strategy for handling these backgrounds with only slight modifications. The Z 
boson background is irreducible, whereas the tibb is reducible and can be reduced so that 
it becomes sub-dominant with respect to the Z boson background. The tibb background 
might seem different from the signal process in terms of its partonic final state, but it can 
mimic our signal, and thus become a relevant background, by “losing” some of the final 
state partons in the detectors. In order to match with the signal’s detector signature, the 
two W bosons originating from the decay of the two top quarks must go undetected, which 
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makes those W bosons the main source of pT in the ttbb background. Although the rate of 
the detector missing the W bosons is expected to be small, the sizable production rate of 
pp —>• ttbb can compensate for this, thus making the pp —)• ttbb process a possibly important 
background. 

A W boson will go unseen in the detector for primarily two reasons: 1) when its 
decay products are not within the experimental acceptance region of the detector due to 
having insufficient pT, supernumerary rj, or both, and 2 ) when its decay products are not 
adequately isolated from other particles, i.e., they are merged with other particles in the 
reconstruction of a given event. For the first case, we define as missed any object that 
satisfies the following criteria: 

• for jets, pT,j < 30 GeV or \r]j\ > 5, 

• for leptons, pT,i < 10 GeV or \rji\ >3 with I = e, p, r. 

In the second case, the following rules determine when a particle is missed: 

• for merging jets, < 0.4 with ji and j 2 denoting any jet pairs including 6 -jets, 

• for merging leptons, ARji < 0.3 with j and I denoting a jet and a lepton, respectively. 

With the acceptance and isolation requirements listed above, we observe that most of the 
background events are from either the fully leptonic or the semi-leptonic decay channels 
of top quark pairs because these channels require that fewer erstwhile visible partons be 
missed in comparison to the fully hadronic top decay channel. 

To devise an adequate strategy for rejecting a large number of these background events 
while preferentially keeping the signal events, we must adopt event selection criteria that 
incorporate the kinematic differences between signal and background event. We observe 
first that, thanks to the heaviness of the parent particles, the signal events will be composed 
of jets that typically have a larger transverse momentum than those found in background 
events. This is a strong hint as to which cuts will be suitable in rejecting the background. 
However, one must be especially cautious in selecting these cuts because the method pro¬ 
posed here is based on extracting which relies in part on the shape analysis of the 
energy distributions for each invariant mass slice. Therefore, cuts should be chosen such 
that they do not considerably distort the energy distributions. For this reason, we pre¬ 
fer using softer cuts than in most searches performed at the LHC, and we choose as our 
baseline selection criteria 


PT,b > 30 GeV, 17/61 <5, ARbb > 0.4, (4.2) 

for identifying the bottom jets in all events that we analyze. In addition, as a good 
approximation, a flat selection rate of 66% is assumed for any single bottom jet to account 
for the 6-tagging efficiency.^® 

^°Note that in practice, this efficiency will have a dependence on pt of the 6-jet, thus modifying the shape 
(i.e., not just the normalization) of the energy spectra. Nonetheless, knowing this functional form of the 
6-tagging efficiency, we can still recover the “true” shape of the 6-jet energy spectra. 
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In order to further suppress the backgrounds from Standard Model processes, we con¬ 
sider requiring that events have a large missing transverse momentum. In signal events, 
the missing transverse momentum is expected to be determined by some combination of 
the new particle masses, and thus will be large. On the other hand, the missing transverse 
momentum in background events is determined by the larger of the total hardness of the 
event, the mass of the Z boson, or the mass of the top quark. Therefore, a large "pT cut 
allows us to efficiently discriminate the signal events from the background. However, in 
our case special care is needed in deciding the scale of this cut, as there is the risk of 
this cut introducing unwanted bias in the energy distributions. In particular, the missing 
transverse momentum can be interpreted as the recoil of the invisible particles against the 
visible, which seems to imply that a large fT cut is likely to select only events with very 
hard visible particles and correspondingly induce some bias in the 6-jet energy spectrum 
toward higher energies. As mentioned before, this could lead to a misidentification the 
value of and as a consequence, an innacurate measurement of the associated masses. 
Fortunately, the relatively large mass hierarchy between the gluino and the neutralino in 
our signal process ensures multiple hard 6-jets on average and thus a sizable recoil for the 
invisible neutralinos. We therefore anticipate that the Ebt, distribution will only be mildly 
affected, even with a fairly hard pT cut. For our signal and backgrounds, we impose 

pT > 200 GeV, (4.3) 

which strongly suppresses the backgrounds with negligible deformation of the Ef,f) distri¬ 
butions. 

In addition to the pT cut, we introduce another cut that requires each 6-jet pT to have 
some minimum angular separation from the pT vector. This enables us to avoid events 
where the measured missing energy is caused by the mismeasurement of jets. For our 
analysis we require 


APifr, I^t) > 0.2, (4.4) 

which has negligible effect on the shape of the Ef^h distributions. 

In Table 1 we show the cross sections for both signal and background events after 
applying the set of cuts listed above. We clearly see that the ttbb background is sub¬ 
dominant with respect to the Z + 6666 background. We also remark that the expected 
signal-to-background ratio (S/B) is large, which is certainly favorable for extracting E'^^ 
from the Ej,f, distribution. Indeed, if new physics particles are discovered in the forthcoming 
runs of the LHC, it would then be natural to discuss measuring their masses in the channels 
where there is a clean signal, and hence a large S/B. In this sense, the context in which 
we present our mass measurement technique is expected to be typical when attempting a 
mass measurement beyond the precision of the order of magnitude, as we do here. 

Other than the above-mentioned backgrounds, QCD multi-jet production pp —)• bbbb is 
another possible source of background events from the SM, in which the missing transverse 
momentum typically arises from imperfectly measuring the energy of jets. Unsurprisingly, 
an accurate estimation of this background is quite challenging because it involves detector 
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effects. We expect that a great deal of the QCD multi-jet background would largely be 
suppressed by the cuts eqs. (4.2), (4.3) and (4.4) to the point that it becomes sub-dominant, 
and in the following we do not taking this background into account (see, for example, 
Ref. [20]). 

4.2.2 Combinatorial background and mixed event subtraction 

As mentioned earlier, the procedure of phase space slicing inevitably requires the formation 
of the invariant mass of two visible objects. Since the process eq. (4.1) under consideration 
has pair-produced parent particles which both decay to indistinguishable visible children, 
grouping the four 6-jets into two pairs gives the correct choice in only one case out of 
the three possible combinations. Following the strategy outlined in Section 3, we form all 
possible pairs and obtain an inclusive energy distribution. We then subtract the contri¬ 
butions originating from the wrong combinations by estimating the corresponding distri¬ 
bution through the event mixing technique as described before. In order to validate the 
performance of this mixed event subtraction scheme in our example, we first study a large 
number of pure signal events where no selection cuts are imposed and without including 
backgrounds. We then discuss how the inclusion of backgrounds complicates the estimate 
of the combinatorial background when using the mixed event subtraction. 

i) Pure signal: For our study it is necessary to check that 6-jet pairs’ energy and invariant 
mass distributions are well reproduced by the mixed event subtraction scheme. In order to 
apply it as per eq. (3.4), we need to obtain the distributions of observables given by forming 
pairs of 6-jets belonging to different events as explained in Section 3. In principle, there are 
several options for choosing the two events that one can use to compute these observables. 
For example, one can compute the observables using all possible pairs of events, meaning 
that each event is reused many times, or one can use a procedure such that each event is 
used only a few times. The detailed way in which the event mixing is done can, in principle, 
affect the result. However, in most cases, the alternative methods give rise to only minor 
differences in the relevant result Among those possibilities, we show results for which 
the events were mixed as follows: given a sample of N events, 

(1) we first randomly shuffle and reorder the events to remove any potential correlation 
between events arising from the way in which the events were generated, 

(2) compute the different-event observables by taking 6-jets in the ith and (i -|- I)th 
events,so that each event is used twice 

(3) finally renormalize these distributions to weigh as much as the contribution from the 
incorrect pairings in the signal sample that we intend to remove, i.e., two thirds of 
the total number of events in the signal sample. 

^^One can also form “events” out of randomly selected sets of four particles from the entire event sample 
and compute the observables by forming pairs from the particles now constituting these new “events.” We 
evaluated the efficacy of constructing the distribution to be subtracted using this alternative method and 
found little difference between the end results, both in the distributions and in the energy and mass values 
extracted. 

^^The last event is mixed with the first one. 
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Figure 2. The left panel shows the di-5-jet invariant mass distributions that are normalized to 
an integrated luminosity of 3 ab~^. The right panel shows the Ri distribution over rribb- 

We label the inclusive invariant mass distribution formed from all the pairs in the 
same event as FsEi'mbb) and we denote as FEEirribb) the distribution obtained from pairs 
in different events. Then, given an invariant mass value, we take from the same-event 
sample all the pairs whose invariant mass lies within the range of interest and plot the 
spectrum of the energy of the sum of the two 5-jets, E^b- We repeat the same operation 
in the different-event sample, and we obtain the E^b spectrum for the hxed ranges of rribb- 
Denoting the spectrum from the same-event pairs as fsEiEbb), the spectrum from different- 
event pairs as foEiEbb), and the resultant spectrum from the mixed event subtraction as 
fs{Ebb), our estimate of the energy distribution from the correct pairs is; 

fs{Ebb) = fsE{Ebb) - foEiEbb) (4.5) 

from which we will ultimately extract the rest-frame energy value E^f^ for each fixed rribb- 
Similarly to the subtraction for the Ebb spectrum, we obtain an estimate of the overall 
invariant mass distribution for the correct pairings using the mixed event technique, which 
is 


Fsirribb) = FsEirribb) - EoEirribb) - 


(4.6) 


Note that this distribution is not used to extract the masses, but serves as an example 
of the effectiveness of the mixed event subtraction scheme. To quantify the fidelity of the 
distribution obtained from the mixed event subtraction scheme, we define a ratio 


Ri = 


Ni,s 

Ni,c 


(4.7) 


where is the i-th bin-count of the “subtracted” distribution Es{mbb) and W,c is the 
corresponding bin-count in the distribution obtained considering only the correct pairs. 
The latter is obtained by exploiting the fact that in simulated events all of the history of 
the particles is available. The left panel in Figure 2 compares the invariant mass distribu¬ 
tion from the correct pairs, shown as the blue dot-dashed histogram, and the distribution 
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Figure 3. The two plots in the left panel show the di-&-jet energy distributions for 300 GeV (top) 
and 700 GeV (bottom) nominal mass slices. They are normalized to an integrated luminosity of 
3/ab. The color codes in the plots in the left panel are the same as those in Figure 2. The two 
plots in the right panel show the respective R distributions over Ebi,. For computing (R), only the 
data within the two black vertical dashed lines is taken into account. 


obtained by mixed event subtraction, shown as the red solid histogram. To show how much 
the original invariant mass distribution is contaminated by the combinatorial background, 
the distribution before the subtraction procedure is also plotted as the green dashed his¬ 
togram. Each bin count is normalized to an integrated luminosity of 3 ab“^. We observe 
that the distribution obtained from the mixed event subtraction is very close to that ob¬ 
tained from the correct pairs. A more quantitative comparison is also provided in the right 
panel of Figure 2, showing the bin-by-bin ratio Ri for the distribution of the events against 
rrihb- All of the bin counts are quite close to their associated theoretical values. In fact, 
~ 1 in all the range of rribb- We do not show the ratio Ri in the vicinity of both kine¬ 
matic endpoints because the significantly smaller Aj,c at the endpoints leads to unreliable 
Ri values. Besides visualizing the effectiveness of the mixed event subtraction, this check 
enables us to confirm that for a given invariant mass slice, a similar amount of data remain 
available after the mixed event subtraction compared to what would have been available 
were we able to completely eliminate the combinatorial background. 
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In Figure 3 we compare the energy distribution from the correct pairs and that ob¬ 
tained from the mixed event subtraction. The two distributions are shown for a mass slice 
275 GeV < rrn,i) < 325 GeV in the upper panel of the figure and for another mass slice 
675 GeV < < 725 GeV in the bottom panel. The ratio of the correct pairs and the 

subtracted histograms is provided for each choice of rrihh as well. Since the fit to extract 
from the energy spectra will be performed using only the data around the peak, we show 
the bin-by-bin ratio only for the energy range that corresponds to the full width at half 
maximum (FWHM), which is indicated by black dashed lines in each plot. In this case, we 
see that the energy spectrum processed with the mixed event subtraction (blue histogram) 
is also quite close to the associated theory expectation (red histogram). To be more quan¬ 
titative, we compute the average of Ri in the FWHM range. This average is denoted as 
{R) and is close to 1, which suggests that the mixed event subtraction scheme works quite 
well, i.e., the shape of the energy spectrum is reasonably preserved. From this exploratory 
analysis, we expect that the extraction of from the energy distribution obtained by the 
mixed event subtraction is unlikely to have major bias due to the subtraction. 

ii) Background and “signal-background interference”: Once the SM backgrounds 
come into play, there is a non-trivial complication that is introduced by the event mixing. 
Since we are not aware a priori whether a given event is from the signal or the background, 
it is not possible to perform the event mixing using only the signal events. Therefore, 
the distribution returned by the whole operation of the mixed event subtraction scheme 
contains “signal-background interference”, i.e., picking one particle from a signal event and 
the other from a background event. In principle, the overall kinematic characteristics of the 
background events differ from those of the signal events, and therefore these interference 
pairs will make the overall distribution deviate from that of the pure signal or pure back¬ 
ground combinatorially-generated background.As a consequence, a naively subtracted 
distribution would be distorted with respect to the distribution of a pure signal sample. 

To understand the quantitative impact of the inclusion of physical backgrounds, we 
first need to assess the hierarchy of the effects that arise from the simple addition of these 
backgrounds and from the event mixing. We focus on the situations where riev events have 
been collected after the application of selection cuts, e.g., eqs. (4.2)-(4.4). These events 
come both from the signal process and from background processes. In general, we have rig 
signal events and = Uev ~ background events.However, in a situation in which a 
mass measurement is attempted, we expect that the signal will dominate the backgrounds. 
Us 3> ni). Under this assumption, we can quantify how likely the event mixing procedure is 
to form pairs where both particles come from the signal process, both particles come from 
the background processes, or one particle comes from signal and the other from background. 
Clearly, the probabilities to select a particle from a signal or a background event in the 

^®For the dominant background in our case (i.e., Z + the pure background combinatorial distribution 
is somewhat tricky; the distinction between correct and incorrect pairings is meaningless because the asso¬ 
ciated event topology is ill-defined. However, in the interest of generality, we imagine that our background 
can also give rise to fictitious correct and wrong combinations. 

^^Since different types of backgrounds, in principle, will form different distributions, we here assume only 
a single type of background to avoid any potential complication. 
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sample are 


Tig Til) 

Ps = -^-, Pb = -^, 

Us + Uh Us Us + rib ns 

for a signal and a backgronnd event, respectively. Therefore, most of the pairs formed in 
the event mixing procednre are made with two particles from the signal process. Pairs 
made of two particles from the background are much less abundant, and in fact arise only 
in a small fraction, n^/n^, of the cases. Strikingly, pairs made with one particle coming 
from the background and one particle from the signal are much more abundant than pure 
background pairs, as their probability is 2 x n^/ns. 

The effect of the pairs involving both backgrounds and signal is not predictable unless 
one specifies the signal. However, some general features of this “interference” contribution 
to the event mixing estimate of the distribution for the correct signal pairings can be 
easily guessed. First, the “interference” distribution tends to produce an underestimation 
of the bin-counts in the estimate eq. (3.4) of the distribution consisting of pairs of 6-jets 
originating from the same decay. The reason is that in the inclusive distribution from all 
pairs in the same event, the first term in eq. (3.4), there are contributions stemming from 
pairs of events coming both from signal or both from background, but there is no way to 
construct a hybrid of the two: 



da 

dEhb 


(all-pairs) 


das das 

dEbb dEbb ’ 


(4.9) 


where by as and as we mean the signal and background contributions, respectively. On 
the other hand, in mixed events, we have 


da 

dEbb 


(different-events pairs) 


dass da SB dasB 

dEbb dEbb dEbb 


(4.10) 


As suggested by Figure 3, the contribution from pairs where both events come from the 
signal, does a good job of estimating the effect of the pairs of 6-jets not coming from 

the the same decay in the signal. Similarly, the contribution from pairs of events where both 
events come from the background, is a good estimate of the combinatorial background 

generated from the background itself. Therefore, the contribution is the piece that 
typically ruins the result because it gives rise to an excessive subtraction in eq. (3.4). 
Obviously, this phenomenon cannot be avoided since we are unable to distinguish signal 
and background events with absolute certainty. The presence of this type of “interference” 
background is inherent to the event mixing technique and dealing with it requires special 
care. A more quantitative argument about the interference is available in App. A for more 
interested readers. 

Given its peculiar origin, it may be desirable in some cases to remove the “signal- 
background interference” contribution. In order to do so, one must discuss the shape of 
this distribution, which in general depends on the signal and therefore is a priori unknown. 
However, some general features of the “signal-background interference” distribution can be 
predicted using the following argument. The distribution that arises from pairs made of 
one background and one signal particle feels in part the kinematics of the signal events and 
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Figure 4. The di-6-jet energy distributions of the true background and the interference for 
300 GeV (left panel) and 700 GeV (right panel) nominal mass slices. They are normalized to an 
integrated luminosity of 3 ab“^. The black distribution is obtained by subtracting the blue one by 
the red one. The vertical black dashed lines denote the associated fitting range for each slice. The 
bottom panel shows the performance of the proposed fitting template for the effective backgrounds. 


in part that of the background events. The background is typically expected to have softer 
particles than the signal, and therefore, the “interference” b-^et pair energy distribution is 
expected to be skewed towards energies that are somewhat larger than the characteristic 
values of the distribution for pairs made out of two background particles. For our example 
process, we display the distribution from pairs of signal particles in Figure 3, while the 
distribution from pairs of background particles and the interference distribution are shown 
in Figure 4. Comparison of these distributions confirms our intuition from the argument 
above. 

From Figure 4 we see that the interference effect dominates the pure background effect 
essentially everywhere in the distribution. As discussed above, this is due to the fact that 
we envision a situation where signal events are much more abundant than background 
events. The shape of the interference contribution is also quite different from that of the 
pure background contribution. 

When discussing results in a later section, we will study the effect of the removal of 
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Figure 5. Average in the FWHM range of the bin-by-bin ratio of the energy distributions from 
the event mixing and from just the correct pairs of 6-jets, (i?) = 1 implies a good match. 

background contributions to our results. With this goal in mind, we study what functional 
form describes the total effect of backgrounds; that is, the distribution from the pure back¬ 
ground pairs plus that from the “signal-background interference” pairs. The combination is 
shown in Figure 4; in the lower panel, we show a possible fit of this distribution. Due to the 
importance of the signal in determining the shape of the “signal-background interference” 
distribution, we decided to model the total effect of backgrounds with a function of the 
family eq. (2.8). The fit result in Figure 4 is rather good, but we do not attach any special 
significance to this finding. In fact, a better description for this background may exist and 
might be preferred. More generally, we stress that the “signal-background interference” 
distribution is not universal, and our choice could be unreliable for other signals. In our 
application to the gluino decay process, the fairly good description provided by eq. ( 2 . 8 ) 
and shown in Figure 4 is satisfactory for our current purposes. 

5 Mass measurement results and discussion 

In this section, we demonstrate the application of the proposed technique to the gluino 
decay. Results on the mass measurement from fitting the energy spectra for the compound 
system of two 6 -jets are presented in the following subsections along with the possible issues 
and limitations of our method. In the final subsection, we discuss possible improvements 
of the mass measurement with the aid of the di-jet invariant mass endpoint. 

5.1 Measurement of gluino and neutralino masses 

Following the strategy outlined in the previous sections, we present results for the determi¬ 
nation of the masses of the gluino and the neutralino from the 6 -jet energies. The energy 
spectra that we study are obtained from simulated event samples generated as described in 
Sec. 4.1 for both signal and dominant background processes at the 14 TeV LHC. We also 
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Figure 6. Comparison of the goodness of the fit between massive template eq. (2.8) (blue) and 
the massless template eq. (2.7) (red) as fitter for the energy spectrum of the 6 -jet pair system for 
various values of the invariant mass of 6 -jet system fhbi). 


recall that the relevant channel is characterized by a large missing transverse momentum 
and four bottom-tagged jets which are selected as per eqs. (4.2)-(4.4). Since the primary 
interest of this paper is to study the theoretical aspects of energy peaks in a multi-body 
decay, rather than data analysis under realistic statistics, we take a sufficiently large num¬ 
ber of events to minimize potential statistical fluctuation within the data sample, which is 
then normalized to an integrated luminosity of 3 ab“^. 

Note that we study the distribution of the sum of the energy of 6 -jet pairs (say, 6 i and 
62 ), i.e., Ehb = Eb^ +Eb 2 , for which the associated invariant mass values belong to a narrow 
range such that 


mbh e [rribb - ^mbbl2, rubb + Ambb/2] , (5.1) 

where Anibb denotes the width of the mass window. We henceforth identify every individual 
mass window by its central value rhbb- Due to the indistinguishable nature of the final 
state particles, we form the Ebb distribution from all possible pairs of 6 -jets in an event, 
and subsequently apply the mixed event subtraction technique described in Sec. 4.2.2 to 
eliminate the contamination from the pairs of 6 -jets not coming from the same gluino. In 
Figure 5, we present the average of the bin-by-bin ratio of the distributions from only 
correct pairs and from the mixed event subtraction technique at various values of fhbb- We 
remark that for each point in the figure, the average is limited to the energy range defined 
by the full width at half maximum (FWHM) of the distribution. The figure suggests that 
the average deviation between the two distributions is, at most, about 8 %, meaning that the 
mixed event subtraction scheme reproduces the original distribution fairly well. Although 
we do not show it here, we also remark that for a given fhbb the standard deviation in the 
bin-by-bin ratio is small enough that the distribution from the mixed event subtraction 
consistently tracks the corresponding distribution from the correct pairs. 
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For each rhbb, the rest-frame energy of the 6-jet pair system (i.e., is extracted from 
the energy distribution by fitting the data to a template function. We have two possible 
template functions, given in eqs. (2.7) and (2.8), and we use both of them to see which one 
better suits the data and to see if there are significant differences between the rest-frame 
energy and E* value found by the two functions. Since eq. (2.7) is suitable only for the 
case where the relevant ifibb is effectively negligible, we expect an increasing discrepancy 
between the results obtained by eqs. (2.7) and (2.8) as fhbb grows. The reduced values 
for each fit to the energy spectrum are shown in Figure 6. This is a measure of how 
well the template function describes the data globally, but we do not necessarily attach 
any statistical meaning to it. We instead use it as a measure for the distance between the 
two template functions. Looking at the hgure, we observe that for all fhbb, the massive 
template describes the data as well as or better than the massless template. As expected, 
the performance of the massless template becomes progressively worse as fhbb increases. 
The massless template also seems inferior to the massive template from another aspect as 
it typically returns an estimate of the rest-frame energy that is larger than the expected 
value. On the other hand, the massive template does not introduce such a pronounced 
bias. 

To demonstrate the difference between fits made with the two templates, we show in 
Figure 7 sample fit results for two different nominal fhbb values, 250 GeV and 650 GeV. In 
the left panels, we provide results from applying the massive template, while in the right 
panels we provide results from applying the massless template. The errors quoted for the 
extracted were estimated at 95% conhdence interval (C.I.) from the variation of the 
of the fit. For both of the fhbb values, we see that the massless template estimates E^^ as 
being slightly larger than the estimate given by the massive template, and the discrepancy 
is larger as we go to larger fhbb- Although the discrepancy is within the 95% confidence 
interval, we feel that this is an important characteristic of the fit results. In fact, the 
massless template has a systematic tendency to return larger E^j^ in general, which implies 
the introduction of a possible bias to the mass measurement. The results of the hts for 
all of the values of fhbb are reported in Table 2, from which we see that massless template 
consistently overshoots the estimate of obtained from the massive template, and that 
it also overestimates the theory values of E^f^. 

Finally, we take all the values of E^^ obtained from the fitting the energy spectrum 
for each fhbb and fit them to the line given by eq. (1.3), taking into account and displaying 
the associated errors from the fit procedure that we used to extract the values of The 
expression in eq. (1.3) can be adjusted to match our example process as follows: 


77'* _ 

^bb — 


m| - + mgfc 

2m?, 


(5.2) 


We perform the fit of eq. (5.2) on both of the results from the massive and massless 
templates. For the fits using the massive template, we use only the results obtained for 
fhbb in the range from 200 GeV to 650 GeV, in which the errors from the ht of the energy 
spectra are quite small. Other choices of the fhbb range give similar mass measurements, but 
we simply make a conservative choice for the range of fhbb included so to avoid the values 
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Figure 7. Sample fit results for extracting using the massless template (right panels) and 
the massive template (left panels). The chosen invariant mass slices are rubb G [225, 275] GeV (top 
panels) and rubb € [625, 675] GeV (bottom panels). We report only statistical errors. Each fit 
range is chosen such that it roughly corresponds to the relevant FWHM. 


of fhbb where fewer events are expected. From the results extracted using the massless 
template, we choose to ht eq. (5.2) only for fhbb in the range from 200 GeV to 600 GeV, 
where it is more reasonably accurate to treat the 6-jet pair system as massless. The fit 
parameters are the slope of eq. (5.2) and its vertical intercept (that is of rhbb = 0) and 
we denote them as 


s = 


1 


2mg 


and 


2 2 
— m% 

y= 3 X 


For our spectrum, the theory values are 

s = 4.2 X 10"^ GeV"\ 


2mg 


y = 595 GeV 


(5.3) 


(5.4) 


Fitting the line eq. (5.2) on the results obtained from the energy spectra, we obtain the 
best-fit lines shown in Figure 8, which correspond to 

s = (4.8 ± 0.3) X 10"^ GeV■^ y = 597 ±5 GeV, (5.5) 

for the massive template and 

s = (5.2 ± 0.3) X lO”'^ GeV“\ y = 606 ±6 GeV, (5.6) 
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rubb 

Fit range 

Theory 

Massive [y^/d.o.f] 

Massless [x2/d.o.f] 

200 

[400, 

1000] 

612.5 

614.8+2l:o [0-064] 

619.1121:2 [0-054] 

250 

[400, 

1000] 

621.9 

627.2t{ll [0.039] 

634.8l}|^ [0.035] 

300 

[400, 

1000] 

633.3 

640.91}^;^ [0.089] 

654.3l}|^ [0.057] 

350 

[440, 

1000] 

646.9 

659.2l}|2 [0.074] 

673.ll}|^ [0-063] 

400 

[440, 

1040] 

662.5 

670.7l}|^ [0.074] 

693.7im [0.061] 

450 

[500, 

1040] 

680.2 

694.81}^^ [0.041] 

715.9l}f^ [0.037] 

500 

[540, 

1040] 

700.0 

716.31}^^ [0.033] 

738.811^1 [0-050] 

550 

[600, 

1100] 

721.9 

742.ll^®;[ [0.038] 

760.0111^ [0-064] 

600 

[640, 

1100] 

745.8 

768.91^^1 [0.026] 

787.9im [0.074] 

650 

[700, 

1200] 

771.9 

802.2+21-0 [0.030] 

810.0ll^;[j [0.063] 

700 

[740, 

1240] 

800.0 

832.712+4^ [0.011] 

840.9lll;o [0.060] 

750 

[800, 

1300] 

830.2 

871.4l?8^^^ [0.017] 

865.31^1? [0.060] 

800 

[840, 

1340] 

862.5 

910.5l??5*g [0.024] 

908.3l^lj [0.067] 

850 

[880, 

1340] 

896.9 

952.81^^2^9 [0.019] 

961.113+^ [0.12] 

900 

[920, 

1400] 

933.3 

998.01^1^ [0.040] 

1015.2l3i;6 [0.21] 


Table 2. The fit results for fifteen invariant mass slices. For each fit, a mass slice of 50 GeV was 
chosen, for example, for fhbb = 200, Ebb is selected such that the corresponding nibb is between 175 
and 225 GeV. The bin size for all energy distributions is 20 GeV. The error estimation for each fit 
parameter is performed by 95% confidence interval. All values but given in GeV. 


for the massless template. Not surprisingly, the values extracted using the massive template 
are closer to the input values than those from the massless template, although even the 
massless template was able to get a rough estimate of the rest frame energy values of 
the bb system for the chosen values of fhbb- In Figure 9, we also show the 68% and 95% 
confidence level contours obtained from the variation of the fit of the slope and intercept 
parameters; the result with the massive template are in the left panel and that with the 
massless template in the right panel. One can clearly see that the distance between the 
theory values and the best-fit values for the case of the massive template (left panel) is 
smaller than that for the case of the massless template (right panel). 

As mentioned before, s and y can be easily converted into the masses of gluino and 
neutralino. Based on eqs. (5.5) and (5.6), we obtain the following measurements of the two 
masses: 


Massive template : rrig = 1042 ± 65 GeV, = —159000 ± 59000 GeV^ , (5.7) 

Massless template : rrig = 964 ± 56 GeV, = —240000 ± 41000 GeV^ . (5.8) 

We remark that the gluino mass, while quite precisely determined, is underestimated by 
about 20%, with the value from the massive template being closer to the true value than 
the value from the massless template. The neutralino mass is poorly determined using 
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Figure 8. The fit of the data points with eq. (5.2). The theoretical expectation for 

the given mass spectrum is represented by the solid black line. The data points obtained by fitting 
the Ehb distributions with massless and massive templates are marked by “cross” and “square” 
symbols, respectively. The mass measurement done with cross symbols is represented by the red 
dashed line. For the blue dot-dashed line, the measurement is done for the data points for which 
the massless template work reasonably well. 


both the massive and the massless template. Possible causes of this poor estimation will 
be discussed in the next subsection. 

5.2 Study of systematic effects 

The results in the previous subsection on the measurement of the gluino and neutralino 
masses are fairly good, considering the challenging circumstances of the mass measurement, 
particularly the fully indistinguishable character of final state particles in our chosen signal 
process. Despite an adequate result for the gluino mass measurement, the neutralino mass 
measurement is very poor; the only conclusion that one is able to draw is that the neutralino 
in our example process is consistent with being massless. 

As noted already, the measurement of for each fhhb is statistically compatible with 
the theory value, but still results in a mass measurement that is systematically overesti¬ 
mated. From the fit of the data in Figure 8, the mismeasurement of primarily implies 
a slope larger than that predicted by theory, which consequently implies that the extracted 
gluino mass is biased towards values smaller than the true mass. This bias is not par¬ 
ticularly worrisome per se, as it is about 10%. However, given the relation between the 
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Figure 9. Contour plots in the plane of s (= l/2mg) vs. y (= (m~ — m|)/2mg) around the best 
fit values for the fit results with massive (left panel) and massless (right panel) templates. 



slope 

steeper 

consistent 

shallow 

y-intercept 

larger 

^g,ext < ^(7,in 

m-o , m?o ■ 

XT,ext Xi,in 

~ ^g,in 
m-o . < m?o ■ 

XT,ext Xi,in 

^ 3 ,ext ^ 3 ,in 

2 2 

mto . ~ rnto ■ 

xV,ext Xi.in 

consistent 

^g,ext in 

m?o , < m-o. 
Xi.ext x¥,in 

ext ^ '^g,m 

9 2 

mlo . ~ ^-0 • 

xV,ext Xi^n 

^ 3 ,ext > ^ 3 ,in 

m?o . > "l-o . 
xV^ext xV.in 

smaller 

^g,ext < ^g,in 

2 2 

m -o . Ri m -o . 
xV,ext XT-™ 

'l^g,ext ~ ^ 3 ,in 

m-o . > m-o . 
xV.ext Xi:in 

^ 3 ,ext > ^ 3 ,in 

m?o . "l-o • 

xV.ext Xi>in 


Table 3. Comparisons of extracted mass parameters with corresponding input values for the nine 
possible combinations of over-, under- or consistent estimation of the slope and the intercept of the 
straight line eq. (5.2) fitted on the {nibb, E^b) data. The orange table cell corresponds to the result 
of the fit of the data in Sec. 5.1. 


masses and the observables in eq. (5.3), it turns out that this underestimation of the gluino 
mass severely affects the neutralino mass determination. More generally, there are nine 
possible cases based on under-, over-, or consistent estimations of the slope and intercept 
of the straight line in eq. (5.2). The implication of each case in terms of the extracted mass 
parameters is summarized in Table 3. 

It is interesting to examine possible causes of this bias in the best-fit line of Figure 8, 
which also serves as a basis for possible improvements of our method. In order to clarify 
the origin of the incorrect estimation of E*, we study the following potential sources of 
inaccuracy in our fits of the energy spectra: 

i) an imperfect fit of the data with the massive template eq. (2.8); 

ii) contamination due to the background; 
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Arubb 

Event mixing 

Background 

Background fit 

Cuts 

OS 

50 GeV 

Yes 

Included 

No 

Yes 

CS I 

50 GeV 

Yes 

Included 

Yes 

Yes 

CS II 

50 GeV 

No 

Included 

No 

Yes 

CS III 

50 GeV 

No 

Not Included 

- 

Yes 

CS IV 

2 GeV 

No 

Not Included 

- 

Yes 

CS V 

2 GeV 

No 

Not Included 

- 

No 


Table 4. Description of the original sample (OS) and several selected check samples (CS). The 
width of the ranges of rrihb for the discretization of the multi-body phase space is reported in the 
first column. Samples marked as event mixed are those in which the mixed event subtraction has 
been carried out. In those marked as “no”, the correct pairs are identified in the event record 
and so eliminate the effect of combinatorial backgrounds. Samples where the background has been 
completely neglected are marked in the third column. For the samples where the background has 
been added, we report in the fourth column if we have added a template to fit the background 
events to the overall fit of the data . Finally, in the fifth and final column we report if selection 
cuts eqs. (4.2) through (4.4) have been applied to the events or not. 


hi) biases introduced by the event mixing subtraction; 

iv) finite size of the rribb range used to discretize the multi-body phase-space; 

v) biases due to events selection. 

For the first potential source, we recall the discussion in Section 2 where the massive 
template function was introduced; this template had a maximum at only when w —)> oo, 
which corresponds to producing the gluino(s) at rest. For practical cases, w is finite and 
the maximum of the function appears at a somewhat larger value than E^b- On the other 
hand, physical energy distributions can have the maximum at in particular, for cases 
where rribb can be treated as effectively massless. Therefore, the relevant fit could result 
in a value that does not match the corresponding expectation. Fortunately for the case at 
hand, we find that w is large enough to cause only a negligible shift in the peak position, 
i.e., such a potential mismatch is very tiny. Consequently, we do not ascribe the systematic 
overestimate E^^ to the inaccuracy of the template function eq. (2.8). 

In order to see the effect of the other four potential sources of bias on the final result, 
we conduct a dedicated analysis for each. In each analysis, we repeat the same procedure 
as described in the previous section, that is to say we extract the values of E^/^ from an 
event sample that incorporates the effect under study. The event samples for the study of 
these possible effects are denoted as “Check Sample” (CS). We then compare the results 
obtained from those Check Samples to those obtained from the Original Sample (OS). The 
attributes of the check samples that we have considered are summarized in Table 4 and 
are also described in the following. 

The first check sample enables us to find the effect of the background on the extraction 
of E^^. We study first the pure background energy distributions in order to calibrate the 
template function describing them in the fit. This calibration is done for each of the rhbb 
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Figure 10. Comparisons of fit results from the five samples used to assess the effects of several 
potential sources of bias in the gluino mass determination as described in Table 4. 

slices. Although there are two types of backgrounds, the dominant SM background and the 
interference from the event mixing, we employ a single template in eq. (2.8) to describe both 
of them collectively. We then repeat the fit of the energy spectra for each fhhb, including 
the template function for the backgrounds as well. The results in the determination of 
for each rhbb are labeled as “CS I” an plotted as red open circles in Figure 10. In this figure, 
the left panel shows the absolute shift of the measured from the corresponding theory 
value for each sample. The right panel shows the ratio of the values of E^^ from the check 
samples and the corresponding value in the original sample. We observe that the effect of 
background modeling is negligible for all fhbb, and from the right panel of Figure 10, we 
can in fact see that this engenders less than 1% of the shift in E^f^. 

The next potential source of bias that we study is the event mixing, which is studied 
by the second and third check samples, denoted by CS II and CS III in the following. In 
these samples we use the event record to identify the correct pairs of jets coming from 
the same gluino, and therefore we obtain the correct energy spectra without applying the 
mixed event subtraction. The two samples CS II and CS III differ by the inclusion of the 
SM background. The results for the determination of are reported in Figure 10 by 
blue filled triangles and blue open triangles, respectively. From the figure we see that the 
determination of is significantly improved. In fact, in the left panel of Figure 10 we 
can see that a mild positive shift of the determined E^^ values still exists, but is greatly 
reduced compared to what we had with the original sample. We remark further that only 
minor differences are found between the results obtained from the check samples CS II and 
CS III, which can be taken as another way of confirming that the effect from background 
events is negligible. Therefore, we conclude that the effect of event mixing is a major cause 
of the shift that we observe in the gluino mass determination. 

Next, we study the effect of the discretization of the nibb spectrum by taking smaller 
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ranges for the fhbb window. This is performed on the check sample denoted as CS IV. 
This narrow rubb range analysis is intended to provide an improvement to the previous two 
analyses where the combinatorial issues were artificially resolved using the information in 
the event record. In fact, the check sample CS IV for this analysis is similar to CS III 
except for the Anibb- The results of this analysis are reported in Figure 10 by black filled 
rhombuses, which suggests that the effect of discretization of the rribb spectrum is negligible. 
This is not surprising in light of the following observation: one can easily figure out that, for 
a given nominal value of fhbb the the from eq. (5.2) for {fhbb + 25) GeV {fhbb — 25) GeV 
is at most about 15 GeV larger (smaller) than the for the nominal fhbb- The absolute 
size of this shift is already quite small and is further reduced by the fact that for each nibb 
range we observe the sum of all the contributions below and above fhbb- In ah; we expect 
a very small net effect, making the discretization of the rubb spectrum in increments of 
50 GeV suitable for the precision sought. 

Finally, we study the bias induced by the selection cuts. To assess their effect, we 
produce a sample along the line of CS IV, but being fully inclusive in the signal phase- 
space. The result of fits performed on the energy spectra from this sample are reported in 
Figure 10 by black open rhombuses. The use of a fully inclusive sample gives from the 
fits that agree with the theory predictions within few percents up to fhbb = 650 GeV. 

For the check samples II, III, IV and V we remark that the agreement of the fit results 
with the theory value deteriorates as one gets closer to the endpoint of the range covered 
by rubb- We observe that good agreement is retained up to fhbb = 650 GeV, which is in 
the falling tail of the m-bb distribution, as apparent from Figure 2. We suspect that the 
mismatch of the fitted and the theory values is connected to the massive template 
becoming less accurate in fitting to the data. Indeed, in Figure 8 we can see that the error 
estimation on the htted becomes larger for fhbb > 650 GeV. In a realistic application 
of our mass measurement method, we would not know up to what precise value of rubb the 
massive template can be trusted. However, it is clear that the values of m-bb for which the 
extracted F^^ from the ht comes with a large error should be avoided. We remark that in 
Figure 8, all the fit results for nibb > 650 GeV have a significantly large error, thus clearly 
signaling a transition to a region of nibb where the fit template eq. (2.8) can no longer be 
trusted. A more detailed investigation of this transition boundary is beyond the scope of 
our paper and we instead refer to Ref. [30] for a more systematic study of it. 


5.3 Improving the mass measurement using the mbb endpoint 

In this section, we discuss a possible improvement of the mass measurement with the aid of 
the kinematic endpoint of the dijet invariant mass distribution. The point is that without 
prior knowledge of the masses, it is not possible to say whether the measurement eq. (5.7) 
has a bias. However, we can devise a check and an improvement of the obtained measure¬ 
ment using an independent observable. To this end, we study a possible combination of the 
result of fits to energy spectra with the measurement of the endpoint of the invariant mass 
distribution of the pairs of 6-jet, which we denote as This observable is a simple 
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Figure 11. The functional dependence of eq. (5.10) for various gluino masses. The rribb endpoint 
is set to be 1.1 TeV. The black solid line denotes the case where rrig is identical to that of our study 
point. 


function of the gluino and neutralino masses: 


^bb =mg-m^, 


(5.9) 


and is expected to be very useful in combination with the results of fits to energy spectra. 
In fact, eq. (5.2), in light of eq. (5.9), can be rewritten as 


TP* — 

^bb - ^bb - 


(m; 


max'\2 -2 
’ ’ I ! f L' ' 


bb ! 


bb 


2m?, 


(5.10) 


which is a straight line in the plane {Elb^ffi^b) described by just one free parameter. This 
should be compared to the previous equation we used to find the masses, eq. (5.2), where 
there are two independent parameters: the slope and the constant term of the straight line. 

If the mhb endpoint is assumed to be well-measured, we can use eq. (5.10) to more ac¬ 
curately fit our results to a line in the plane (m^^, This relation is shown in Figure 11 
for = 1.1 TeV. The determination of kinematic endpoints is a well-established tech¬ 

nique, as experimental collaborations have adopted the idea for precision measurement (e.g. 
Ref. [9] ). Several phenomenological works have developed techniques of detecting kinematic 
endpoints or edges in one-dimensional data [43, 44] and even in multi-dimensional data [45]. 
To construct the template for the endpoint extraction, we closely look at the dijet invariant 
mass distribution. As explicitly shown in Figure 2, the combinatorial background is a major 
challenge, and we again conduct the mixed event subtraction with a signal-background- 
combined event sample. The left panel of Figure 12 exhibits the distribution constructed 
by all possible pairs (green dashed histogram) and the distribution after the mixed event 
subtraction (red solid histogram). To guide the location of the theoretic endpoint, we lay 
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Figure 12. The left panel shows the di-6-jet invariant mass distributions that are normalized to 
an integrated luminosity of 3 ab“^. Both signal and background events are plotted with selection 
cuts in Section 4.2.1 imposed. The right panel shows the endpoint extraction by fitting the rubb 
data between 1000 and 1200 GeV with the template in eq. (5.11). 

a black, vertical dashed line. Since the full shape of the mbb distribution in this three-body 
decay is well-known (see, for example. Refs. [46, 47]), one could employ the relevant ana¬ 
lytic expression to fit the signal portion. However, considering potential shape distortion 
by the mixed event subtraction, non-trivial spin effects, matrix element etc. (although not 
significant), we instead focus on local information in the vicinity of the kinematic endpoint 
and describe it by a straight line. We also observe that the residual distribution beyond 
the endpoint can be also accommodated by a single straight line. We therefore introduce 
the following template g{mbb) to extract the maximum 

g{rnbb) = si{mbb - to^“) + ^ 2^66 + c, (5.11) 

where c, and are fit parameters responsible for slopes, y-intercept, and the 

endpoint, respectively. We remark that similar strategies have been implemented by ex¬ 
perimentalists (e.g.. Ref. [9]) and in many phenomenological studies (e.g.. Ref. [4]). 

Our fit result with the above template is reported in the right panel of Figure 12. The 
data points (blue dots) in-between 1000 and 1200 GeV are taken for the fit, and the best-fit 
lines are represented by two red straight lines. The extracted endpoint value is 

= 1112.1 ± 3.5 GeV, (5.12) 

where the quoted error is again a statistical uncertainty. The result is quite precise as 
expected. Using this extracted and the relation in eq. (5.10) as a template for 

the fit of the data points in Figure 8 along with eq. (5.9), we obtain a mass 

measurements: 

rUg = 1231 ± 30 GeV, = 119 ± 30 GeV. (5.13) 

In this case, as well as for the analysis in the previous sections, the fit is performed between 
Tfibb = 200 and fhbb = 650 GeV, with the data points obtained from the massive template. 
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This result is more accurate and in better agreement with the expected values than what 
we obtained in eq. (5.7) using only energy distributions. Therefore, we conclude that, 
depending on the accuracy with which the rrihh endpoint can be experimentally determined, 
the addition of information from the rribi, endpoint can grant a very significant improvement 
to our results obtained from only the energy spectra. 

6 Summary and Conclusions 

In this paper, we have discussed how to use the energy spectra of visible decay products 
for the measurement of masses of “parent” particles in semi-invisible multi (i.e., more than 
2)-body decays. The results are an extension of previous results regarding the properties of 
the energy distribution of the massless decay products in a two-body decay. In particular, 
we extended the results for two-body decays by discretizing the multi-body phase-space 
and considering it as the combination of multiple two-body decays, as suggested by the 
recursive factorization formula for the multi-body phase space. 

The fictitious two-body systems that are involved in the recursive factorization of 
phase-space are necessarily massive. Therefore, we utilized the results of a companion 
paper [30] on the description of energy spectra of massive decay products in two-body 
decays. Armed with these results, we should be able to fit the energy spectra of the visible 
part of the fictitious two-body decay and extract an estimate of the masses involved in the 
process using the results of these fits. 

A particularly challenging aspect of our analysis had to do with the fact that the 
decaying particle (the particle of interest) is typically produced in association with other 
particles. In this situation, it is possible that some of the “child” particles from the decay 
of the parent are identical to those contained in the rest of the event, which means that 
in the process of reducing the multi-body final state of the parent decay to one with fewer 
bodies, the particle pairings that we perform may unintentionally include particles which 
have nothing to do with the mass measurement at hand. In particular, the parent particles 
are often produced in pairs: if the two parents in each event undergo the same decay 
process, then it is clear that this combinatorial background is inevitable. These particles 
extraneous to the decay potentially hamper the mass measurement of the parent, and thus 
the contamination they add must be addressed. Most of the general discussion above can 
be succinctly illustrated by the consideration of a suitable example. Furthermore, tackling 
a concrete example enables us to quantify the quality of the mass measurement that one can 
achieve using our method. With these goals in mind, we studied in detail the production 
of a pair of gluinos in a supersymmetric model, where R-parity is conserved, and in which 
the gluinos directly decay to bbxi (a three-body decay) via an ojf-shell bottom squark. To 
this end, we simulated events for the process 

PP ^ 99 ^ bbbb + pT ■, 

including the relevant SM contribution. We identified selection cuts to remove the SM 
backgrounds to a level that further clears a path toward the successful determination of the 
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masses of the new physics states. Simultaneously, we attempted to minimize the changes 
in the shape of the energy spectra caused by these event selection criteria. 

Our actual analysis then starts with forming all possible pairs of bottom quarks in each 
event that passed the selection in eqs. (4.2)-(4.4), from which we then obtained a distri¬ 
bution of the energy of the b quark pairs, Ebf,. This distribution includes the contribution 
from pairs formed by bottom quarks not originating from the same decay, i.e., the wrong 
combinations. Note that the final state of each decay in our chosen process is made of two 
(visible) indistinguishable particles, and so the pollution from the combination of particles 
not coming from the same decay is even more severe; in particular, each event gives 6 
combinations of two bottom quarks, out of which only 2 are correct, cf. the case of distinct 
particles a and b from each decay, which would give 2 correct ones out of 4 combinations. 
To remove this adulteration of the event sample, we subtracted an estimate that we ob¬ 
tained using the event mixing technique described in Sec. 3. This estimate was obtained 
from pairs of b quarks taken from different events, and we showed that the contribution 
from pairs of b quarks not coming from the same gluino can thus be effectively removed, 
as seen in Figure 3. In other words, this method has a natural tendency to maintain the 
shapes of the distributions that we need to analyze to carry out our measurement. 

The energy spectra, once effectively rid of the contribution from pairs of b quarks 
coming from different gluinos, were fitted in a region around their peak with the function 
eq. (2.8), which is taken from our work Ref. [30] and briefly described in Sec. 2. For 
comparison, these energy spectra were also fitted with the function eq. (2.7), which was used 
in our original paper on the peak of energy spectra of massless particles. The comparison 
of the results from fits with the two functions highlights the improvement achieved by our 
new result for massive decay products. The better description of the energy spectra with 
the function eq. (2.8) can be seen in the comparison of the for the various fits performed, 
as reported in Figure 6. 

The result of the fits of the energy spectra is the extraction of the function parameter 
which is exactly the energy of the system of the two b quarks in the rest frame of 
the gluino that generated the two b quarks. The determination of Ef^^ is the core of our 
procedure as this value is connected to the masses of the gluino and the neutralino via 
eq. (5.2): 


- m2 + rhl^ 

- 2mg ' 

for a pair of b quarks of mass rhbb- The results of the extraction of Ef^ for several choices 
of the mass of the two b quark system were shown in Figure 8. The determination of Eff^ 
for each rubb was fitted using the straight line eq. (5.2) given above. This fit is essentially 
our mass measurement, as the gluino mass corresponds to the slope of the line and the 
neutralino mass to the constant term of the straight line. The resulting mass measurement 
was given in eq. (5.7), which was found to be within 20% of the gluino mass and a rather 
poor determination of the neutralino mass. This inaccuracy of the mass measurement is 
mainly due to the fact that in each fit of energy spectra, we tend to overestimate the energy 
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In Sec. 5.2, we studied several possible causes for this error and (in the end) identified 
the modest shape changes due to the mixed event subtraction as the primary source. 

In order to improve the mass measurement, we then studied how including information 
about the endpoint location in the mbb distribution altered the quality of the mass measure¬ 
ments. If well-measured, this quantity should correspond to the mass difference between 
the gluino and the neutralino, and hence is expected to aid in the measurement of these 
two masses. Indeed, a much more accurate measurement was obtained using information 
from this observable, as shown in eq. (5.13). 

As the analysis we have discussed in this work used parton-level events, we finally make 
a brief comment on the potential impact of the effects that would be necessary to include 
in a fully realistic study beyond this level [e.g. parton showering, hadronization, detector 
response, etc.). Since our target final state is the bottom quark jet, present experience with 
such objects in the LHC experiments make us confident that our conclusions would not 
change significantly if we went beyond a parton-level analysis. In fact, we can look at similar 
analyses with 6-jet energy spectra and notice that detector effects and the improvement of 
event generation beyond the parton-level induce small effects that are well under control 
for the level of accuracy of the measurements discussed in this work. For example. Ref. [19] 
performed a detector-level study on the 6-jet energy spectrum from the top quark decay to 
measure the top quark mass, and found that the extracted E* is in good agreement with the 
corresponding input value within the associated uncertainty. More dedicated study on the 
higher order effects in the 6-jet energy spectrum from top quark decay has been conducted 
in Ref. [48] and demonstrate that these effects are negligible for our purpose here. Finally, 
the CMS collaboration adopted the proposal in Ref. [19] and measured the top quark mass 
using the 6-jet energy spectrum [49]. The measured top quark mass from the extracted E* 
was quite consistent with the world-average value; in particular, despite the uncertainties 
associated with 6-jet energy, a percent-level precision measurement of the top quark mass 
was obtained. From all these considerations, we anticipate that for the signal discussed 
in this work the reconstructed 6-jets do not develop a substantial systematic bias in their 
energy distribution, and the relevant uncertainties can be under control. 

In conclusion, we have demonstrated the use of energy distributions, and in particular, 
the region close to their peaks, to measure the masses of new physics particles involved 
in a single-step multi-body decay. Taking the example of gluino production and decay 
in a R-parity conserving supersymmetric model, we have found that using only visible 
decay products of the gluino decay g —)• 66%, it is possible to measure with good accuracy 
the masses of both the gluino and neutralino, with the best results being obtained when 
both the energy and the invariant mass distributions of the pairs of 6 quarks are used. 
Rather strikingly, the mass measurement technique that we discussed does not actually 
use information about the missing momentum, except only for event selection purposes. 
The example that we considered also required a proper removal of the effect from pairs 
of 6 quarks not coming form the same gluino. To this end, we have shown that the event 
mixing technique is especially well-suited. We anticipate that our general methodology 
here can be applied to mass measurements in other processes. 
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A Event mixing techniqne with a signal plus background sample 

We begin with the total number of pairings of 6-jets (denoted by N]\fc) and with the data 
sample consisting of Ug signal events and ni, background events. Remembering the fact 
that there are six possible pairings out of four 6-jets in each event, we have 

Nisfc = 2ns + 2nb + Aus + 4nb (A.l) 

where the last two terms represent the total number of wrong pairings of 6-jets denoted by 

Nwc = 4(ns -b Ub) ■ 


As mentioned before, our unawareness of which are signal and which are background 
events precludes us from having the event mixing (symbolized by (8>) only between signal 
events or background events. Therefore, if we perform the event mixing procedure on 
the total [us -b Ub) events, we then have signal-signal mixing, signal-background mixing, 
and background-background mixing. Since there are four 6-jets in each event, a single 
event mixing enables us to have up to 4 = 16 6-jet pairings. Of course, for practical 
purposes, one could use a subset of these 16 pairs, say m pairs. Now that m is given by a 
common prefactor for every single event mixing, we then need to know the number of event 
mixings. The total number of signal-signal mixing is evaluated by the two-combination in 
the binomial coefficients: 


Us ® Us 



and so is that of background-background mixing; 


nb®nb = 


UbjUb - 1 ) 

2 


(A.2) 


(A.3) 


Likewise, the number of signal-background mixing, i.e., interference, is expressed as follows: 


Us ( 8 ) Ub 



UgUb- 


(A.4) 
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Suppose that we use m mixed pairings out of 16 mixed pairings in each event mixing. 
Denoting by Nmc the sum of eqs. (A.2), (A.3), and (A.4), we eventually use m- Nmc ^-jet 
pairs to estimate the distribution from the N\yc wrong pairs formed in the same event 
distribution.^^ Since in general N]yc / each mixed pair should be re-weighted in 

making the “different events” distributions, so to match the contribution of Nwc wrong 
pairs. Keeping this issue of the normalization in mind, let us first see how the wrong 6-jet 
pairings in the signal events can be treated by the mixed event subtraction even in presence 
of a small background. Since the total number of mixed pairings is normalized to Ny/c^ the 
presence of the background affects the weight of the “different event” distributions made 
by just signal events. In fact, in the “different event” distribution, once normalized so as 
to match the contribution of the Nyyc wrong pairs, the fraction of 6-jet pairs stemming 
from the signal-signal mixing is given by eq. (A.2) times a rescaling factor to match the 
normalization to Nyyc: that is 


risins^ _ 4(n. + m) __ A _ (A^) 

2 n,(n,-l)/2 + n,nb + nbK-l)/2 ^ 

where the common prefactor m is omitted for simplicity and the approximation is done 
with the assumptions of Ug ^ and ;:?> 1. The implication of this result is that the 
combinatorial background (caused by the signal itself) can be almost completely eliminated 
with the event mixing technique thanks to the dominance of the signal assumed throughout 
in this paper. In the total “different event” distribution reweighed as to match the expected 
number of wrong pairing Nwc we find that the fraction of 6-jet pairs from the background- 
background mixing is 


UbiUb - 1) 


4{ns + Ub) 


Using - l)/2 UgUb + Ubiub - l)/2 


Xug Ug 


(A.6) 


In the same distribution we find that the fraction of 6-jet pairs from the signal-background 
mixing is 


UgUb X 


A{ng -t- Ub) 


Using - I)/2 UgUb -k Ubiub - l)/2 


8nb ( 1 - — -k — 
Ug Ug 


(A.7) 


Comparing eq. (A.6) and eq. (A.7), we conclude that the effect of “interference” is, in 
general, more important than the contribution from the background-background mixing, 
as also argued in the main text. Besides corroborating the argument given in the main 
text, these equations quantify more precisely the effect of background, which becomes 
more important when the mass measurement strategy described in this paper is applied to 
situations where S/B \s less favorable than that in our example process. 


^®Note that the only thing that we know is -|- nj,, neither Us nor rib separately. 
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